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I. INTRODUCTION 

When colloidal particles are solvated in a polar fluid they usually release counterions 
and therefore they acquire an electric charge. Then, a solution of heavy charges (the 
colloidal particles) and light charges (the counterions) is obtained. A model of such 
charged colloidal suspensions in sedimentation equilibrium based on the local density 
approximation has been proposed^ and numerical computations about that model show 
that the heavy charges accumulate at the bottom of the container while the light 
particles accumulate at the top. The suspension behaves like a condenser: a vertical 
electric field is induced into the suspension. Also, several experimental measurements 
of the densities profiles in such suspensions have been done. 

In this paper we consider a simple model of sedimentation equilibrium, taking into 
account only the Coulomb interactions between the particles and the gravitational 
force (in particular we neglect excluded volume interactions), which is exactly solvable. 
The model is a two-dimensional two-component plasma in a gravitational field, that 
is a system composed of two species of particles with electric charges itq and masses 
M±. The Coulomb interaction is logarithmic in two-dimensions, and the particles are 
submitted to a uniform vertical gravitational field. 

The two-component plasma model has been solved by Gaudin^ when the inverse 
temperature (3 verifies F := = 2, and much work has been done with that model. 
Studing the electrical double layer for the two-component plasma, Cornu and Jancovici^ 
introduced a very general method for treating the two-component plasma in presence 
of an external potential. We use here their method to deal with the gravitational field. 

In the following section we present the model and recall some results of the method 
as given in Ref. ^ In section III, we solve the model and compute quantities such as 



the density profiles and the electric field in the suspension. 



II. THE MODEL AND THE GENERAL METHOD OF RESOLUTION 

The model is a two-dimensional system of particles of charges ±5 and masses M±. 
The position r of a particle is represented by its Cartesian coordinates {x,y). The 
gravitational field is g = —gy. The particles are in a container of height h and infinite 
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width. The bottom of the container is at y = 0. The half-spaces y < and y > h are 
impenetrable to the particles. The interaction of the particles with the gravitational 
field is 

M+gY,yt + M,gJ2yr (2.1) 

i j 

where yf^ is the altitude of the i^^ particle of charge itg. Expression ( p.ip can be 
rewritten as 

Meff 9 (y: vt + E 2/7 j - 9^efr {y. yt - E yj^ (2-2) 

where Mes = (M+ + M_)/2 and qEes = -{M+- M_)g/2. This shows that the system 
is equivalent to a system where all the particles have the same mass Mefj but there 
is an external electric field E'eff acting on the particles. For this reason it is useful to 
define k± = f3M±g, ko = (A;+ + A:_)/2 and 6k = — k-. 

For F > 2 a point-particle model is unstable against the collapse of pairs of particles 
of opposite charge. To prevent this collapse we introduce some short range cutoff by 
representing the particles as hard discs of diameter R and obtain results in the small- 
R limit. In fact, at F = 2, only the one-body density diverges, while the n-body 
correlations (n > 2) have well-defined limits as ^ 0. From now on we shall consider 
that the temperature is such that F = 2. The model is solved in the grand canonical 
ensemble. It is convenient to introduce position-dependent rescaled fugacities that have 
inverse length dimensions: m±{T) = mQexp{—k±y). There are two length scales in this 
problem: the gravitational length given by k^^ and the screening length which in fact 
is giveni by the inverse of the rescaled fugacity ttiq^. Here we concentrate on the usual 
case close to the continuous medium limit where uiq ^ ko. 

Using the equivalence between the two-dimensional Coulomb gas at F = 2 and the 
free Dirac fieldiii, one can write the grand partition function in the form 



det 



/ N 1 + , ,1 — (Jz 

+ m+lr) h m_(r) 



(2.3) 



2 ' ' 2 

where ^ = axdx + (^ydy is the Dirac operator in two dimensions and ax, cry and cr^ 
are the Pauli matrices. The particle densities ps (where s = ±1) and the truncated 

(2)T 

two-body densities p),J, can be obtained from the Green functions Gss' defined by 



^ + m+(r) — - — + m_(r) — — 



G = 16{r - r') . (2.4) 



2 



One finds 



Psir) = ms(r)Gss(r,r) , 



(2.5) 



and 



p2r(r,r') = -ms{r)mAr')GsAr,r')Gs's{r',r). (2.6) 

It is convenient to write the fugacities as m±(r) = m(r) exp[— (ib)2y(r)] where 
m(r) = moexp[— A:oy] and ^(r) = dty/A, and to define 

gss'{r,r') = exp[-sF(r)]G,,. (r, r') exp[-sV(r')] . (2.7) 

Then, in terms of the operators a = dx + idy + dxV + idyV and = —dx + idy + dxV — 
idyV, the functions Qss' are given by the partial differential equations 

m(r) + {m{r))^^ a g^^{r, r') = 6{r — r') , 

m(r) + a {m{r))^^ g (r, r') = 6{r — r') , 



and 



f-+(r,r') = - [m(r)] ^ a(/++(r, r') , 
5+-(r, r') = [m(r)]"^ g — (r, r') . 



(2.8) 
(2.9) 



(2.10) 
(2.11) 



III. CALCULATIONS AND RESULTS 



A. The Green functions 

Because of the translational invariance along the x axis we can suppose that x' = 0. It 
is useful to work with the Fourier transform ^^s' of Green functions gggi defined by 



1 /■+°° 

gss'{r,r') = — / g.^'iv^v' ,k)e"''' dk . 



(3.1) 



From equation ( |2.8| ) we find that the Fourier transform g^^ satisfies the ordinary 
differential equation 

\k + 6k/i)ko -{k + 6k/ if - mle-^'^'^y] g++{y, y', k) 

+kodyg++{y, y', k) + d^g++{y, y', k) = -moe-''oy6{y - y') . 

In terms of u = exp[— A^oy], u' = exp[— fcoy'] and fss'{u,u',k) = gss'iy,y' ik), we have 
for /_|_-|_ the equation 



(3.2) 



dlf++{u,u',k) 



f++{u,u,k) 



mo 



6{u - u) , (3.3) 



where v\ = [k — {k- /2)Y /k^. A similar equation is obtained for / by replacing v'^ 

by = [-k - {k+/2)f/kl. Equation Q can be solved in terms of the modified 



Bessel functions and Ky^^: 



J++[u,u ,k) = —Vuu' 
ko 



. ko ko ko ko 



(3.4) 

where u< = min(n, n'), u> = niax(u,u'), and A and B are constants (with respect to 
u) that depend on the boundary conditions. Note that the sign of can be arbitrarly 
chosen; we choose positive. Using equation ( 2.10| ) we find 



/ u , kj = —I- — vuu 

ko 



AU^+^C^u) - BK,^+,C^u) (3.5) 
ko ko 



+ 



-/,^(^n')i^.++.(f^n), ifn'<n. 
where a is the sign of A; — (fc_/2). 

For a source point in the allowed region < y' < /i, in the impenetrable region 
y < the particle fugacities vanish and the equations for the Fourier transforms of 
G++ and G ^ reduce to 

{k + dy)G++{y,y\k) = Q, ik-dy)G-+{y,y',k) = 0, (3.6) 
with the solution 

G++(y,y^A;) = Ce-^^ G.+ iy,y' ,k) = De'^v . (3.7) 

Similar equations hold for y > h. 

The boundary conditions are that ggg' must vanish when y — > ±00 and be continuous 
at y = and at y = h. The boundary conditions when y — > — cxd and equations ( p.7| ) 

lead to /++(!, u',fc) = if /c > and / \-{l,u',k) = if A; < 0, while the boundary 

conditions when y — > +00 give at y = h: f++{uh,u' , A;) = if A; < and / \^{uh, u' ,k) = 

if A; > 0, where we have defined Uh = exp(—koh). Using these boundary conditions 



in (3.4) and ( p.5| ) we find the values of A and B which depend on the sign of k and the 
sign of k — k-/2. We shall use the following notation: A^, (u', z/_|_) and i?^, (n',z^+) are 
the values of A and B when k has sign a' and k — k^/2 has sign a. Using this notation 



we have 



,rnoM^u')K,^„{^u^,)+K,{^u')h^^{^u^,) 



A^l„.l ,.\— jr^ ^"'U ^ '^^ fco ^ '^"T" ^ fcO ' ----v feO " ^--^-T" V ftp (r,0„\ 



A'L{u',v) = K,+,{-^ 



(3.8b) 
(3.8c) 



Similar calculations can be done for / Then changing the integral over k in (| 

into an integral over u we find for the Green functions 



(3.8d) 



gss{r,r') = ^ exY>[isk^sX /2]\/ uu' 

r"+00 



2/ Iy{^^u<)Ky(^u>) cos{kQXu) dv 
Jo ko ko 



+ 



+ 



+ 



+00 







ko ko 



,mo 



iskoxv^^ (3.9) 



ko ko 

A-_{u\ u)hC^u) + BZiu', u)K,C^u) 
ko ko . 



-iskoxu 



dv 



lk-s/2ko 

while g |- and are given by dl^) and (|ll|). 

An interesting limit is when the height h of the container goes to infinity. One can 
take the limits of A'^,{u' , u) and B^,{u' , v) when u/j — > in ( [3.91 ). One finds 



5fss(r,r') = ^exp[is/c_sx/2]\/'u.'u' 



+00 



,mo 



Iu{- — u<)Ky{- — n>) cos{koxu) dv 

ko 



ko ko 



+ 



+ 



r lX^u)IA^u')'^f^e-">'«-'i. (3.10) 

Jk^j2ko ko ko ^-i(fc^) 



mo . ,mo .mo. 
-U{—u)Iy{—u )Ky{—) 
ko ko ko 



+ < 



Iu{-;—u)Ky{—u) + Ky{—u)Iu{—u) Ky{ — ) 
ko ko ko ko J ko 



+Ky{—u)K^{—u )Iu{-rr) ) 

ko ko feo vr 



^-"^ ko 



-dv 
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B. The charge density 



The densities ps are given by equation ( ^.51) . The first integral in ( ^.91) and in (|3.10| ) 
diverges for u' = u and must be cut-off as mentioned in section^ an equivalent method 
is to introduce a cutoff t'max in the integral. The other terms are finite when n < 1 
(y 7^ 0). The charge density \s p = q{p+ — p-) and has a finite limit when the cut-off 
R goes to zero. Figures ||, § and ^ show the charge density profile when the heavy 
particles have charge +g, and for mo/kQ = 10, /c_|_ = l.bko, /c_ = 0.5/co and different 
heights h of the container {koh = 0.5, 3, and 10 respectively). 

When h is of the same order or less than the gravitational length k^^ we notice 
a strong accumulation of heavy particles at the bottom of the container and of light 
particles at the top, while in the middle region the charge density is almost zero. In 
that zone there is an almost constant electric field. The suspension behaves like a 
condenser. This was noticed before in other model^. Note that we have considered 
the usual physical case where niQ ^ k^. In the hypothetical case where k^ ^ mo the 
neutral zone would not exist. We shall assume from now on that mo ^ ko. 

When h ^ ko^ , there still is a strong accumulation of heavy particles at the bottom 
of the container but the light particles no longer accumulate on the top of the container. 
There is an excess of light particles at some altitude that will be shown later to be of 
order k^^ ln{mo / ko) . 

Let us consider the case h = +oo. The charge density can be written as 

K 1 ('^) 



p(r) = q 



mlu^ /-fe+Zafeo 



2tt 



fc_/2A:o 



I,A—u] 



lu- 



ko ko 



2Iv[—u)K^{—u)Ky[ — ) 
ko ko ko 



mo , 
ko ' 



2sin7rz^ 



TT 



T (IM. 



(3.11) 



dv . 



First we shall study the charge density near the bottom of the container and at inter- 
mediate altitudes. Remember that we are in the case mo ^ ko- If y ^ k^^ ln(mo/A;o) 
and n > 1/2 then using the asymptotic expansion of the Bessel functions we have 



mot 



47r 



■ exp 



-A;oy- 2^(1-6^^0^) 
ko 



(3.12) 



From this expression we notice that for y <^ kQ^ln^mo/ko), the function q6k p{y) is 
a decreasing function of y. Furthermore from expression (|3.12| ) we see that the layer 
of heavy particles has a thickness of order mg ^, the screening length. Now at high 
altitudes for y ^ ln{mo/ko) we have 



p{y)^-(l^^l -/ shi{'Kv)K^{^uf du . (3.13) 



m-Q-U^ /■fc+/2A:o rk-/2ko 

The asymptotic behavior when y ^ oo of the integral J^"^'^^'^ sm{7ri')Ki^{^uY dv is 
computed in the Appendix. One finds if ^4. 7^ and /c_ 7^ 



2fc>/fco 

p(y) gsgn(*)-^--2- 



mo 



^ ^ ^ ^ ^ ^ 

/cor(A;>/2A;o)^ sin(7rA;>/2A;o) ^ , (3.14) 

y 



where A;< = min(A:+,A;_) and A;> = max (/;;+, A;_). From expression ( p. 14 ), we see that 



q3c p{y) is now an increasing function of y. This shows that at an altitude of order 
k^^ ln(mo/A;o) we a local minimum of the function qSk p{y) 

The electric field in the region where the system is almost neutral {m^^ ^ y ^ 
/Cf^^ ln(mo/fco)) can be estimated as follows. Remember that the system is equivalent 
to a system of particles with the same mass in an external electric field E^s = — (^+ — 
M-)g/2q. The screening properties of this equivalent system will ensure that its total 
electric field is zero, so the electric field in our system, created by the accumulation of 
particles at the bottom and top of the container (or at an altitude of order k^^ ln(mo/A;o) 
if /i > A;^^ ln(mo/A;o)), will be equal to -E^s = (M+ - M^)g/2q. 

C. The density profiles 

In the almost neutral region we can estimate the individual densities of the particles 
through the following macroscopic argument. At altitude y the Coulomb gas pressure 
p compensates the gravitational force: 

^ = -gip+M+ + p_M_). (3.15) 
dy 



This equation, the fact that p+ = p^, and the equation of stateH for the two-dimensional 
two-component plasma for F < 2 

Pp = (p+ + p_) (1 - ^) , (3.16) 
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yield the barometric law 

p+{y) = p-{y) = /5oexp 



-ko 1 



-1 



y 



(3.17) 



This can be verified in the exact expressions of the densities given by the Green func- 
tions. The dominant term in the Green function g^^ in that region is 



5++(r,r') ~ gl^{r,r') 



mo 



vr 



■ exp 



ik^x/2]Vuu' 



+ 00 



iu\-, — u^)Ky[- — n> j cos ukox av , 
ko ko 

(3.18) 



and g is the same as 5++ except for the phase factor exp[iA;_x/2] that is changed into 

exp[— iA;_|_x/2], which is in fact irrelevant for the physical quantities. The integral over v 
can be approximated by a sum over integers by using the Euler-MacLaurin expansion 
when mo ^ ko- Then the sum can be performed by using the addition theorems for 
the Bessel functions. This gives, for the first two terms of the expansion, 



g^^(r,r')~ — ^ exp[ifc_x/2]e" 



ko 



^ko(ix-y) 



-koy' 



(3.19) 



+ l^oi^ exp[-koy])Ko{^ exp[-fcoy']) + 
6 fco fco 



As stated in section ||[ when the diameter R of the particles goes to 0, the density 
diverges, so we must cut-off expression ( 3.19| ) when |r — r'| ^0 replacing |r — r'| by i? 
for |r — r'l < R. Then the regularized form of the particle density is 

m(r)^ 



P+(r) = P-(r) 
If moR <C exp(— fcoy) then 



27r 



-Koim{r)R) 



p+{v) = p^{v) ^ poe-^^^y ^ 



(3.20) 



(3.21) 



where po = mQKo{moR)/{2'iT) is the density when g = 0. Notice that the expres- 
sion ( 3.2C| ) of p± is the same as the one for po except that the fugacity m(r) depends 
on the altitude y. We have found a barometric law ( |3.21 ) with mass M-|_ -|- M_ which 



agrees with ( 3.17 ) when F = 2. 

Far away from the neutral zone, for y ^ kQ^ \n{mo/ko) the dominant term of ps in 
the case mo ^ ko is 



Psiy) 



sm.{TTv)Ky{-—u) dv . 
ko 



(3.22) 



The asymptotic behavior of this integral for y ^ oo is computed in the Appendix. One 
finds if fcs 7^ 



2k-s/ko 

87r2 



mo 



ko 



ks/ko ^ Q-kay 



koTik-s/^koY sin(7rA;__ j2A:o) . (3.23) 

y 



The density ps decays Hke ex.p{—PMsgy). The argument of the exponential is the 
same of the one of an ideal neutral gas, but the decay of the density is faster than the 
one of an ideal gas (because of the factor y~^). 

In real colloidal suspensions the mass of a counterion is much smaller than the mass 
of a colloidal particle. If we assume in our model that the negative couterions have mass 
M_ = then the expression ( |3.1[l| ) for the positive colloidal particles is simplified: the 
last integral in ( 3.10| ) vanishes. The dominant term of the profile density of the colloidal 
particles given by ( |3.20| ) is now valid not only in neutral zone but also up to y ^ cxd. 
For the counterions the profile density in the neutral zone is still dominated by ( 3.20| ) 
but for y ^ kQ^ ln(mo/fco) expression ( |3.23 ) is no longer valid because = 0. Instead 
we have (see the Appendix for details) 

We recover an already known result: in this particular case where k- = 0, in the 
considered region, because of the exponential decay of the density p+{y), one can con- 
sider that the system is composed only of particles with charge —q without mass: the 
system is a one-component plasma without a background. The asymptotic behavior 



given by (3.24) was found for the solvable model of the one-component plasma without 



a background^. The following macroscopic argument also leads to expression (|3.24| ): 
the one-component plasma without a background can be seen as a neutralized one- 
component plasma (where we have added a background charge density qp-) immersed 
in an external charged background —qp~. At altitude y the gradient of the thermal 
pressure of the neutralized plasma compensates the electric force exerced on the exter- 
nal background: 

^ = -qp_{y)E{y) , (3.25) 
dy 

where p{y) is the thermal pressure of the neutralized one-component plasma and E{y) 
is the electric field created by the charge distribution —qp^{y). The Poisson equation 
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IS 

^^ = -27rg/5_(y). (3.26) 

These two equations together with the equation of state for the neutralized one- 
component plasma 

Pp=(i- 5) P- , (3.27) 
yield the differential equation for p_ 

1 dp-{y) 



T\ d 



P-{y) dy 



27rTp_iy), (3.28) 



with the solution 



1 -r/4 1 -r/4 , , 

n (y) = — '—— ~ — '—. (3.29) 



which agrees with (|]2D for r = 2. 



Let us now consider that the mass of the counterions is not zero but that <^ k-i 



In that case equation (3.23) becomes for the counterions 



P~{y)-2 . (3.30) 

47r y 

It is interesting to notice that this expression is the same as the one which would be 
obtained for a one-component plasma without a background but with massive particles. 
For a one-component plasma, it is equivalent to consider that the particles have a mass 
M_ or that they are in an external electric field E = —M^g/q and do not have a mass. 
The calculations done in Ref. ^ can be easily extended to the present case and they 
give the same asympotic behavior ( |3.30| ) for the particle density. 

IV. CONCLUSION 

We have discussed the density profiles in this model of sedimentation equilibrium. Un- 
fortunately our model cannot account for excluded volume interactions which are also 
important in the sedimentation equilibrium. Nevertheless, for heights of the container 
of the same order as the gravitational length, we find the same qualitative results as 
in other models^: heavier particles accumulate at the bottom while lighter particles 
accumulate at the top. For heights of the container much larger than the gravitational 
length, we find an interesting phenomenon: lighter particles do not accumulate any 
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longer at the top of the container but there is an excess of them forming a cloud cen- 
tered at an altitude of order of the gravitational length times some logarithm. An open 
problem is whether such a cloud of light charged particles might be observed, because 
we have not considered possible instabilities in the horizontal direction, with the excess 
charge going to the lateral walls for instance. 
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APPENDIX 

We wish to find the asymptotic behavior of the integral J^"/^'^" sin{TTi')K^{mQu/kQ)'^ dv 
when u — > (y — > oo). We define z = niou/ko, a = ks/2ko, and 

I{z) = r sin{^u)Ky{zf dv . (A.l) 
Jo 

We use the formula 

roo 

Ky{zf = 2 Ko{2z cosht) cosh{2vt) dt , (A.2) 

^0 

in equation ( [A.l| ) and perform the integral over u to find 

^. , /""^ vr — vr cosfavr) cosh(2a;t) + 2t sinfavr) cosh(2at) , , x , ^ 

I(z)=2 ^ — ^ \ \ o ^ — ^ -Ko(2zcoslit)dt. (A.3) 

Jo + 



If a < 1 (ie. if A;_s 7^ 0) then doing the change of variable s = 2zcos\it in (A.3) one 
finds that 

lim z^" ln(l/z)/(z) = 2-^ sin(a^) / s^'^-^Kois) ds . (A.4) 
2^0 Jo 

The integral over s can be done and finally we have 

/(.)^^^-r(a)^sin(avr)^. (A.5) 
This expression is not valid when a = 1, that is in the case when one of the masses of 



the particles is zero. In that case equation (A.3) becomes 



I{z)=4j 2 _^ 4^2 -^o(^^ cosh t) cosh^ t dt . (A.6) 
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Then doing again the change of variable s = 2z cosh t one finds 



I{z) - ^ \ , ■ (A.7) 
^ 2-*o 4 z2ln(l/2;)2 ^ ' 
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FIG 1: The charge density p{y) for mo/ko = 10, k+ = 1.5ko, k- = 0.5ko and 
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FIG 2: The charge density p{y) for mo/ko = 10, k+ = 1.5ko, k- = 0.5ko and h = S/cq 
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